library(ltm) #for cronbach's alpha
library(ggeffects) #graphs
library(ggExtra) #graphs
library(cowplot) #graphs
library(patchwork) #graphs
library(broom) #graphs
library(ggplot2) #plotting
library(dplyr)
library(tidyr)
library(tidyverse)
library(corrplot) #factor analysis
library(psych) #factor analysis
library(factoextra) #factor analysis
library(FactoMineR) #factor analysis
library(doBy) #multicollinearity
library(sjPlot) #plotting
library(viridis) #plotting
library(stargazer) #table output


load("RGC_Time_Experiment_Data.rdata")

#Figure 5####

#ensure analysis only uses people who answered support question
blackorg <- subset(timedv_final, blacksupport >= 0)
nonblackorg <- subset(timedv_final, whitesupport >= 0)

m1giveblack <- lm(QuestionsBlack ~ rgcmean, data = blackorg) #regression for number of
#questions answered for respondents randomized to see the black organization
m1supportblack<-lm(blacksupport ~ rgcmean, data = blackorg) #regression for support of
#organization for respondents randomized to see the black organization

m2givenonblack <-lm(QuestionsWhite ~ rgcmean, data = nonblackorg) #regression for number of
#questions answered for respondents randomized to see the race-neutral organization
m2supportnonblack <-lm(whitesupport ~ rgcmean, data = nonblackorg) #regression for support of
#organization for respondents randomized to see the race-neutral organization

#Table A.5####
stargazer(m2supportnonblack, m1supportblack, m2givenonblack, m1giveblack, type = "html",
          out = "Figure 4 Models.doc")

#get estimate, standard error, p.value, and confidence intervals for regression giving to black organization regression
m1giveblackout <- tidy(m1giveblack, conf.int = TRUE, conf.level = .95)
#remove intercept
m1giveblackout <- subset(m1giveblackout, term == "rgcmean")
#name group 
m1giveblackout$group <- "Predicted Value of RGC on Answered Questions 0-10"

#get estimate, standard error, p.value, and confidence intervals for regression supporting black organization regression
m1supportblackout <- tidy(m1supportblack, conf.int = TRUE, conf.level = .95)
#remove intercept
m1supportblackout <- subset(m1supportblackout, term == "rgcmean")
#name group 
m1supportblackout$group <- "Predicted Value of RGC on 0-10 Support"

#combine two datasets and name them
combined<-rbind(m1giveblackout, m1supportblackout)
combined$estimation <- combined$group
combined$groups[combined$group == "Predicted Value of RGC on 0-10 Support"] <- "  "
combined$groups[combined$group == "Predicted Value of RGC on Answered Questions 0-10"] <- " "
combined$organization <- "Black Organization"

#do the same above process for the nonblack organizations
m2givenonblackout <- tidy(m2givenonblack, conf.int = TRUE, conf.level = .95)
m2givenonblackout <- subset(m2givenonblackout, term == "rgcmean")
m2givenonblackout$group <- "Predicted Value of RGC on Answered Questions 0-10"

m2supportnonblackout <- tidy(m2supportnonblack, conf.int = TRUE, conf.level = .95)
m2supportnonblackout <- subset(m2supportnonblackout, term == "rgcmean")
m2supportnonblackout$group <- "Predicted Value of RGC on 0-10 Support"

combinednonblack <- rbind(m2givenonblackout, m2supportnonblackout)
combinednonblack$estimation <- combinednonblack$group
combinednonblack$groups[combinednonblack$group == 
                          "Predicted Value of RGC on 0-10 Support"] <- "  "
combinednonblack$groups[combinednonblack$group == 
                          "Predicted Value of RGC on Answered Questions 0-10"] <- " "
combinednonblack$organization <- "Race-Neutral Organization"


combinedplot_final<-rbind(combined, combinednonblack)
combinedplot_final$Organization <- combinedplot_final$organization


#begin plotting
p <- ggplot(data = combinedplot_final, aes(x = groups, y = estimate, ymin = conf.low, ymax = conf.high))
p<-p + 
  geom_hline(yintercept = 0, color = "black", linetype='dashed') +
  geom_pointrange(aes(color=Organization),
                  position = position_dodge(width = 0.35), size=0.5)+
  coord_flip() + 
  scale_x_discrete(labels = c("RGC on 0-10 Questions Answered", "RGC on 0-10 Support")) +
  scale_y_continuous(breaks=c(-5,-4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10), 
                     limits = c(-5, 10)) +
  theme_classic() +
  theme(legend.position='bottom', legend.text=element_text(size=12))+   
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(panel.background = element_rect(fill = 'white', color = 'black')) + 
  labs(title = "RGC and Time Experiment", 
       caption = "",  
       x = "", y = "0-10 Support/0-10 Questions Answered") + 
  scale_color_manual(values=c("black", "gray")) +
  guides(color=guide_legend(nrow=2))

p
ggsave(filename="figure5.jpeg", device="jpeg", dpi=500)



#cronbach's alpha of rgc variable####
myvars <- c("ridimp", "selfimage", "impself", 
            "black.fate", "autonomy", "ftblack",
            "close","discrim","rrdiscrim",
            "keptback", "power", "collective",
            "votediff", "blackorg", "black.rights", 
            "workhardindv", "linkfate")

ndexp <- timedv_final[myvars]

cronbach.alpha(ndexp, na.rm = TRUE)

#Table A.15####
#get the means
summary(timedv_final$blacksupport)
summary(timedv_final$whitesupport)
summary(timedv_final$QuestionsBlack)
summary(timedv_final$QuestionsWhite)

#t.test for individual means and confidence intervals
t.test(timedv_final$blacksupport)
t.test(timedv_final$whitesupport)
t.test(timedv_final$QuestionsBlack)
t.test(timedv_final$QuestionsWhite)

#t.tests to see the difference between the two means
t.test(timedv_final$blacksupport, timedv_final$whitesupport)
t.test(timedv_final$QuestionsBlack, timedv_final$QuestionsWhite)

#correlation between supporting the organization and donating to the organization
cor(timedv_final$blacksupport, timedv_final$QuestionsBlack, use="complete.obs")
cor(timedv_final$whitesupport, timedv_final$QuestionsWhite, use="complete.obs")



#Figure A15####
#Lowesss Plots


p_questions <- ggplot(timedv_final, aes(x = rgcmean)) +
  geom_jitter(aes(y = QuestionsBlack, colour = "QuestionsBlack")) + #spreads the points to reduce overplotting
  geom_smooth(aes(y = QuestionsBlack, colour = "QuestionsBlack"), 
              method = loess, se = TRUE, linetype = "solid") + #linear line
  geom_smooth(aes(y = QuestionsBlack, colour = "QuestionsBlack"), 
              method = lm, se = FALSE, linetype = "dashed") + #loess line
  geom_jitter(aes(y = QuestionsWhite, colour = "QuestionsWhite")) +
  geom_smooth(aes(y = QuestionsWhite, colour = "QuestionsWhite"), 
              method = loess, se = TRUE, linetype = "solid") + #linear line
  geom_smooth(aes(y = QuestionsWhite, colour = "QuestionsWhite"), 
              method = lm, se = FALSE, linetype = "dashed") + #loess line
  theme_classic() +
  labs(title = "Feedback", 
       x = "RGC",
       y = "0-10 Questions Answered") +
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_x_continuous(limits = c(0,1)) +
  scale_y_continuous(limits = c(0,10))+
  scale_colour_manual(values = c("red", "blue"),
                      labels=c('Black Organization', 'Race Neutral Organization'),
                      name = "")
  
p_support <- ggplot(timedv_final, aes(x = rgcmean)) +
  geom_jitter(aes(y = blacksupport, colour = "blacksupport")) + #spreads the points to reduce overplotting
  geom_smooth(aes(y = blacksupport, colour = "blacksupport"), 
              method = loess, se = TRUE, linetype = "solid") + #linear line
  geom_smooth(aes(y = blacksupport, colour = "blacksupport"), 
              method = lm, se = FALSE, linetype = "dashed") + #loess line
  geom_jitter(aes(y = whitesupport, colour = "whitesupport")) +
  geom_smooth(aes(y = whitesupport, colour = "whitesupport"), 
              method = loess, se = TRUE, linetype = "solid") + #linear line
  geom_smooth(aes(y = whitesupport, colour = "whitesupport"), 
              method = lm, se = FALSE, linetype = "dashed") + #loess line
  theme_classic() +
  labs(title = "Support", 
       x = "RGC",
       y = "0-10 Support") +
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_x_continuous(limits = c(0,1)) +
  scale_y_continuous(limits = c(0,10)) +
  scale_colour_manual(values = c("red", "blue"),
                      labels=c('Black Organization', 'Race Neutral Organization'),
                      name = "")
library(ggpubr)
figurea10 <- ggarrange(p_support, p_questions, ncol = 2, common.legend = TRUE, legend = "bottom")
figurea10

#Table A.11, A.13 Descriptive statistics of time experiment####
##Table A.11 - descriptive statistics for time study####
median(timedv_final$income, na.rm = TRUE)
median(timedv_final$education, na.rm = TRUE)
mean(timedv_final$female)
median(timedv_final$age, na.rm = TRUE)
mean(timedv_final$democrat, na.rm = TRUE)
median(timedv_final$ideology, na.rm = TRUE)

##Table A.13 - descriptive statistics by condition####
#descriptive statistics for Black organization condition
median(timedv_final[timedv_final$Condition == "Black Organization",]$income, na.rm = TRUE)
median(timedv_final[timedv_final$Condition == "Black Organization",]$education, na.rm = TRUE)
mean(timedv_final[timedv_final$Condition == "Black Organization",]$female)
median(timedv_final[timedv_final$Condition == "Black Organization",]$age, na.rm = TRUE)
mean(timedv_final[timedv_final$Condition == "Black Organization",]$democrat, na.rm = TRUE)
median(timedv_final[timedv_final$Condition == "Black Organization",]$ideology, na.rm = TRUE)

#descriptive statistics for white organization condition
median(timedv_final[timedv_final$Condition == "Race Neutral Organization",]$income, na.rm = TRUE)
median(timedv_final[timedv_final$Condition == "Race Neutral Organization",]$education, na.rm = TRUE)
mean(timedv_final[timedv_final$Condition == "Race Neutral Organization",]$female)
median(timedv_final[timedv_final$Condition == "Race Neutral Organization",]$age, na.rm = TRUE)
mean(timedv_final[timedv_final$Condition == "Race Neutral Organization",]$democrat, na.rm = TRUE)
median(timedv_final[timedv_final$Condition == "Race Neutral Organization",]$ideology, na.rm = TRUE)

#Factor Analysis####
##Figure A.17 - Correlation Plot####
cor.table = cor(ndexp, use = "na.or.complete")
corrplot(cor.table, method = 'shade', order = 'AOE', diag = FALSE)

##Table A.23 - PCA analysis, 3 factors, varimax rotation, standardized loadings####
principals<-principal(na.omit(ndexp))
principals$values

fit <- principal(na.omit(ndexp), nfactors=3, rotate="varimax")
fit # print results

###factor 1: RIDIMP, Selfimage, Impself, ftblack, close
###factor 2: blackfate, keptback, linkfate, blackrights, rrdiscrim, autonomy, discrim (equal btwn 2/3)
###factor 3: blackorg, votediff, collective, power, workhard

res.pca <- PCA(ndexp, graph = TRUE)
#Scree Plot - taken out of appendix
fviz_eig(res.pca, addlabels = TRUE, ylim = c(0, 50)
         )

##Table A.25 - factor regression####
timedv_final$affectidentity <- (rowMeans(timedv_final
                              [c("ridimp", "selfimage", "impself", "ftblack", "close")], 
                              na.rm=TRUE))
timedv_final$commonfatediscrim <- (rowMeans(timedv_final
                              [c("black.fate", "keptback", "linkfate", "black.rights", 
                              "autonomy", "rrdiscrim", "discrim")], na.rm=TRUE))
timedv_final$collectiveact <- (rowMeans(timedv_final
                              [c("blackorg", "votediff", "collective", "power", 
                              "workhardindv")], na.rm=TRUE))
timedv_final$blackorganization <- 0
timedv_final$blackorganization[timedv_final$Condition == "Black Organization"] <-1
mfactor1<-lm(Questions ~ affectidentity*blackorganization, data = timedv_final)
mfactor2<-lm(Questions ~ commonfatediscrim*blackorganization, data = timedv_final)
mfactor3<-lm(Questions ~ collectiveact*blackorganization, data = timedv_final)

summary(mfactor1)
summary(mfactor2)
summary(mfactor3)

stargazer(mfactor1, mfactor2, mfactor3, type = "html", out="factorregresstime.doc")



#Figure A.19/Table A.28 - interaction regression####

#Model

timedv_final$blackorganization[timedv_final$Condition == "Black Organization"] <- 1
timedv_final$blackorganization[timedv_final$Condition == "Race Neutral Organization"] <- 0

#make a scale for income to range from 0 to 1
timedv_final$incomescale <- (timedv_final$income -1) / 15

m_int<-lm(Questions ~ rgcmean*blackorganization + incomescale, data = timedv_final)
summary(m_int)

#Plot
p<-plot_model(m_int, type = "pred", terms = c("rgcmean", "blackorganization"), 
              colors = "bw", ci.style = c("whisker"),
              axis.title = "Questions Answered (0-10)",
              title = "Time Study Interaction Model"
)
p<-p+labs(x = "RGC")+
  scale_y_continuous(limits = c(-1.75, 6))+ 
  theme(
  panel.grid.major.y = element_blank(),
  panel.grid.minor.y = element_blank(),
  panel.grid.major.x = element_blank(), 
  panel.grid.minor.x = element_blank() 
) + 
  theme(legend.position='bottom')+ 
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(panel.background = element_rect(fill = 'white', color = 'black'))

p


#Table A.21 Regression Results with Controls, Meancentered RGC and Reweighted Gender####
timedv_final$rgcmeancentered <-timedv_final$rgcmean - 0.6841852 #the mean of rgc is 0.6841852

timedv_final$regressionweight[timedv_final$female == 1] <- 2.145594  #reweighting for gender, since more women then men in the dataset
timedv_final$regressionweight[timedv_final$female == 0] <- 1

m1meancentered<-lm(Questions ~ rgcmeancentered*blackorganization + 
                     incomescale + female, data = timedv_final)
summary(m1meancentered)

meancenteredreweight<-lm(Questions ~ rgcmeancentered*blackorganization, 
                         data = timedv_final, weights = regressionweight)
summary(meancenteredreweight)

stargazer(m_int, m1meancentered, meancenteredreweight, type = "html",
          out = "Interaction Regression Models and Controls.doc")

#Table A.16 - Multicollinearity####
#need the regression with mean centered rgc, gender, income, experimental condition,
# and experimental condition*rgc mean centered

car::vif(m1meancentered)




